function [ Fgrid ] = Bivariate_4pointRefinement_boundary( grid , omega )
%-----------------------------------------------------------------
% Input:
% grid(m,n) - a grid over Z^2
% Output:
% Fgrid(2*m-1,2*n-1)
%
% ASSUMPTION: m,n > 4
%-----------------------------------------------------------------
% ABSTRACT
% Tensor product, includes boundaries. the middle point is insert in the order y->x of
% calculation
%-----------------------------------------------------------------
% NIr Sharon, 26-05-12
%-----------------------------------------------------------------

m = size(grid,1);
n = size(grid,2);

Fgrid = zeros(2*m-2,2*n-2);


%-----------------------------------------------------------
% First main loops: For vertices (new and old) on the grid

for x=1:(m)
    for y=1:(n)
        % Interpolation
        Fgrid(2*x-1,2*y-1) = grid(x,y);
        % Y direction
        if (y~=n)
        if (y==1)
            Fgrid(2*x-1,2*y) = FourPointAvarage_boundary([grid(x,y) , grid(x,y+1) , grid(x,y+2) , grid(x,y+3)],omega);
        else
            if (y==(n-1))
               Fgrid(2*x-1,2*y) = FourPointAvarage_boundary([grid(x,y+1) , grid(x,y) , grid(x,y-1) , grid(x,y-2)],omega); 
            else
               Fgrid(2*x-1,2*y) = FourPointAvarage([grid(x,y-1) , grid(x,y) , grid(x,y+1) , grid(x,y+2)],omega); 
            end
        end
        end
        % X direction
        if (x~=m)
        if (x==1)
            Fgrid(2*x,2*y-1) = FourPointAvarage_boundary([grid(x,y) , grid(x+1,y) , grid(x+2,y) , grid(x+3,y)],omega);
        else
            if (x==(m-1))
               Fgrid(2*x,2*y-1) = FourPointAvarage_boundary([grid(x+1,y) , grid(x,y) , grid(x-1,y) , grid(x-2,y)],omega); 
            else
               Fgrid(2*x,2*y-1) = FourPointAvarage([grid(x-1,y) , grid(x,y) , grid(x+1,y) , grid(x+2,y)],omega); 
            end
        end
        end
    end
end

%Fgrid(2*m-1,2*y-1) = grid(x,y);

%-----------------------------------------------------------
% Second main loops : grid face centers 

for x=1:(m-1)
    for y=1:(n-1)
        % Y direction        
        if (y==1)
            m_y = FourPointAvarage_boundary([Fgrid(2*x,2*y-1) , Fgrid(2*x,2*y+1) , Fgrid(2*x,2*y+3) , Fgrid(2*x,2*y+5)],omega);
        else
            if (y==(n-1))
                m_y = FourPointAvarage_boundary([Fgrid(2*x,2*y+1) , Fgrid(2*x,2*y-1) , Fgrid(2*x,2*y-3) , Fgrid(2*x,2*y-5)],omega);
            else
                m_y = FourPointAvarage([Fgrid(2*x,2*y-3) , Fgrid(2*x,2*y-1) , Fgrid(2*x,2*y+1) , Fgrid(2*x,2*y+3)],omega);
            end
        end
        % X direction
        if (x==1)
            m_x = FourPointAvarage_boundary([Fgrid(2*x-1,2*y) , Fgrid(2*x+1,2*y) , Fgrid(2*x+3,2*y) , Fgrid(2*x+5,2*y)],omega);
        else
            if (x==(m-1))
                m_x = FourPointAvarage_boundary([Fgrid(2*x+1,2*y) , Fgrid(2*x-1,2*y) , Fgrid(2*x-3,2*y) , Fgrid(2*x-5,2*y)],omega);
            else
                m_x = FourPointAvarage([Fgrid(2*x-3,2*y) , Fgrid(2*x-1,2*y) , Fgrid(2*x+1,2*y) , Fgrid(2*x+3,2*y)],omega);
            end
        end
       % dif(x,y) = m_x-m_y;
        Fgrid(2*x,2*y) = (m_x+m_y)/2;
    end
end
  
end



